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ABSTRACT 


Research  was  conducted  to  determine  the  feasibility  of  adaptively  identifying  the  values  of 
the  hydrodynamic  parameters  of  a  submarine  or  model  from  submerged  free-flight  trajectory  data. 

The  extended  Kalman  filter  approach  which  was  used  was  based  on  the  equation-error  parameter 
identification  principle  instead  of  the  more  conventional  response  error  method,  resulting  in  a  relatively 
modest  computational  requirement.  If  the  controls  are  driven  by  appropriate  pseudo-random  noise 
sequences,  it  is  possible  to  determine  almost  all  the  parameters  to  good  accuracy  from  a  five  minute 
trajectory  sample  which  is  contaminated  by  realistic  additive  noise  plus  bias.  Several  configurations 
of  suitable  instruments  for  measuring  the  trajectory  are  described  and  the  algorithms  for  processing 
the  trajectory  data  are  discussed. 
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I.  INTRODUCTION 


PURPOSE  AND  SCOPE 

The  purpose  and  scope  of  this  research  effort  may  be  briefly  defined: 

Purpose:  Immediate 

To  demonstrate  the  feasibility  of  identifying  the  values  of  the  hydrodynamic 
parameters  of  a  submarine  or  model  from  free-flight  data. 

Ultimate 

To  construct  and  provide  to  the  U.S.  Navy  an  instrumentation  package  and  a  data- 
processing  method  which  are  capable  of  identifying  submarine  or  model  parameters 
in  an  operational  environment. 

Scope:  The  scope  of  the  research  reported  herein  was  limited  to  the  immediate  purpose, 

noted  above.  The  research  consisted,  broadly,  of: 

a.  analytical  study  to  adapt  to  this  purpose  the  parameter-identifying  procedures 
which  we  had  developed, 

b.  to  evaluate  their  performance  by  digital  simulation,  and 

c.  to  indicate  configurations  of  instruments  which  could  provide  the  essential 
input  signals. 

BACKGROUND 

The  system  parameters  identification  problem  has  been  treated  by  many  researchers  over  the 
years  from  various  points  of  view.  The  work  of  Kalman^2  )  has  sparked  a  revolution  in  the  develop¬ 
ment  of  solutions  for  state  estimation  problems.  In  this  report  the  theory  is  extended  to  parameter 
identification  in  a  manner  that  is  inverted  from  standard  applications.  Generally  the  equations  of 
motion  form  the  process  upon  which  measurements  are  made.  The  technique  advanced  here  interprets 
the  equation  of  motion  as  the -measurements,  and  the  error  structure  of  the  true  measurements  as 
the  process;  this  method  is  related  to  the  Equation-Error  Adaptive  Identification  technique. 

One  of  the  more  significant  elements  of  this  program  was,  in  a  sense,  completed  before  the 
study  began.  This  was  the  decision  that  the  identification  method  should  be  based  on  the  equation 
error  rather  than  the  response  error  adaptive  parameter  identifier  technique.  The  motivation  for  this 
choice  is  outlined.  The  equation-error  method  permits  decoupling  of  each  of  the  six  equations  of 
motion  from  the  others;  the  response  error  method,  related  to  conventional  Kalman  filters,  does  not. 
This  fact  enables  the  equation  error  identifier  to  be  far  simpler  in  its  computational  complexity. 

We  briefly  consider  the  relative  complexity  of  the  two  systems.  Neglect  in  each  case  the  problems 
of  noise  and  bias  in  the  instruments,  and  assume  that  all  signals  may  be  measured. 

Each  of  the  six  coupled  nonlinear  equations  of  motion  has  25  (or  more)  unknown  parameters. 
A  response-error  model  must  therefore  be  a  nonlinear  coupled  net  of  dynamical  order  of  six.  In 
addition,  there  are  required  150  nonlinear  coupled  nets  to  calculate  gradients  plus  150  nonlinear 
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coupled  algorithm  nets  to  calculate  the  parameters.  All  these  nets  are  dynamically  coupled  and 
nonlinear,  so  that  the  total  system  is  of  dynamical  order  of  306  (or  more).  The  problems  of  feasibility, 
of  computability,  of  stability,  and  of  uniqueness  of  solution  can  only  be  imagined.  On  the  other 
hand,  if  the  equation  error  method  is  chosen,  each  of  the  six  equations  may  be  treated  separately 
from  the  others;  one  linear  model  without  dynamics  can  be  formed  for  each  equation,  and  the  25 
signals  required  to  estimate  the  gradient  are  available  directly  from  each  error  net.  The  parameter 
estimating  algorithms  require  25  integrators.  Six  separate  linear  problems  of  dynamical  order  of  25 
therefore  result.  Uniqueness  and  the  effects  of  noise  may  be  precalculated. 

In  view  of  the  relative  computational  complexity  of  the  response-error  approach,  and  the 
relatively  advanced  state  of  the  art  of  inertial  sensors,  it  was  concluded  that  this  study  should  be 
restricted  to  the  problem  of  developing  and  modifying  an  equation-error  technique  to  be  applicable 
to  the  problem.  The  result  is  therefore  an  equation-error  Kalman  combination  which  is  an  extension 
of  the  Kalman  principles  in  a  novel  and  unconventional  direction. 

An  appendix  provides  a  detailed  heuristic  explanation  of  the  computational  procedure  for 
those  who  wish  to  use  the  method,  or  to  understand  the  mathematical  techniques.  It  is  not  essential 
for  those  who  wish  to  appreciate  the  results  of  this  research  which  are  presented  in  Sections  VI  through 
XI. 
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II.  CONCLUSIONS 


1 .  It  is  feasible  to  estimate  the  entire  set  of  hydrodynamic  parameters  of  a  submarine  from  free- 
flight  data  with  state-of-the-art  instrumentation  and  an  algorithmic  procedure  which  com¬ 
bines  an  adaptive  equation-error  parameter  identification  procedure  with  Kalman  state 
estimation.  This  approach  results  in  relatively  modest  computational  requirements. 

2.  Errors  of  parameter  identification  in  the  presence  of  instrumentation  noise  range  from  a  few 
percent  for  parameters  which  strongly  influence  the  submarine  trajectory  to  100%  for  those 
which  negligibly  influence  the  trajectory  and  thus  can  be  obscured  by  noise. 

3.  Successful  parameter  identification  requires  excitation  not  only  of  the  bow  and  stern  planes 
and  rudder  but  also  requires  blowing  ballast  and  changing  speed. 

4.  The  state  estimation  and  parameter  identification  algorithms  are  capable  of  estimating  the 
bias  on  the  signals  from  the  various  sensors.  In  addition,  they  can  be  extended  to  enable  es¬ 
timation  of  nonstationary  parameters  and  bias.  These  algorithms  are  thus  capable  of  very 
wide  application  through  the  entire  area  of  the  physical  sciences. 

5.  The  combination  of  parameter  identification  and  state  estimation  techniques  yields  an  extended 
Kalman  filter.  The  problems  of  instability  frequently  observed  in  extended  Kalman  filters 
have  been  fully  resolved  in  this  approach;  it  is  believed  that  this  method  may,  therefore,  have 

a  significant  influence  on  the  practical  application  Of  extended  Kalman  filters  and  on  the  di¬ 
rection  of  their  development. 
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III.  RECOMMENDATIONS 


1 .  Much  more  efficient  computational  procedures  can  and  should  be  developed;  those  used  in 
this  project  were  selected  to  ensure  feasibility,  rather  than  to  achieve  efficiency. 

2.  The  present  parameter  identification  and  state  estimation  algorithms  should  be  extended, 
minimizing  the  uncertainties  of  the  submarine  mass  and  inertias,  the  center-of-gravity  and 
center  of  buoyancy  locations,  to  correcting  errors  due  to  nonconstant  error  sources,  such 
as  sensor  drift  and  Schuler  loop  phenomena,  and  to  considering  the  implications  of  un¬ 
certainty  of  the  ballast  state. 

3.  The  design  of  a  specific  parameter  identification  equipment  should  now  be  considered. 
ORGANIZATION  OF  THE  REPORT 

It  is  assumed  that  the  reader  has  some  familiarity  with  vector-matrix  notation  and  discrete 
Kalman  filtering  theory.  The  report  begins  by  introducing  notations  and  conventions  to  facilitate 
presentation  of  the  basic  ideas  (Section  IV).  An  estimation  problem  for  which  a  solution  can  be 
advanced  is  then  briefly  formulated  in  Section  V.  Section  VI  develops  the  theoretical  bases  for  the 
computer  algorithm  by  first  expressing  the  submarine  estimation  problem  in  the  r.otational  frame¬ 
work  of  Kalman  filtering  theory.  The  general  algorithm  is  then  presented  and  discussed  in  terms  of 
its  stability  and  convergence. 

Section  VII  provides  the  philosophical  motivation  for  submarine  parameter  identification 
simulation  studies;  these  are  summarized  in  Sections  VIII  and  IX. 

Methods  by  which  the  dynamic  trajectory  variables  defining  the  motion  of  the  vehicle  can 
be  instrumented  are  presented  and  discussed  in  Section  X.  The  remaining  sections  (XI  and  XII)  present 
normalized  plots  of  the  computer  algorithm’s  convergence  to  the  solutions  for  the  hydrodynamic 
coefficients. 
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IV.  EQUATIONS  OF  MOTION 


To  facilitate  presentation  of  the  basic  ideas,  vector  notations  will  be  employed  to  express  the 
general  equations  of  a  submarine. 

Detailed  expressions  for  the  force-moment  equations  which  characterize  the  motion  of  a  sub¬ 
marine  are  given  in  *  and  will  not  be  completely  reproduced  here.  It  will  suffice  to  represent  them  in 
the  following  form : 

fi  (xd)  =  ^  r i  (xd>  8 )  +  hj,  i  -  1 , 2,  .  .  .  6 

where , 

i  is  the  equation  number  designating  the  axial,  lateral,  normal,  roll,  pitch,  and  yaw, 

force  and  moment  equations  respectively, 

xd  is  the  dynamic  state  vector  whose  elements  are  the  components  of  vehicle  velocity, 

acceleration,  angular  rates  and  angular  acceleration  in  body  fixed  coordinates 
respectively, 

Cj  is  a  vector  whose  elements  are  the  hydrodynamic  parameters  for  equation  i, 

5  is  a  vector  denoting  the  deflections  of  the  control  surfaces, 

hj  designates  the  thrust  and  weight-blown  effects  plus  other  force-moment  phenomena 

that  act  on  the  submarine, 

r j  (x^,  6 )  denotes  a  function  vector  for  equation  i  whose  elements  couple  the  hydrodynamic 
parameters  into  the  equations  of  motion, 

fj  (xd)  represents  the  mass  or  moment  of  inertia  terms,  accordingly, 

T  superscript  T  implies  “Transpose”. 

To  illustrate  the  notation  a  shorter  form  of  the  axial  force  equation  (i  =  1 )  will  be  written  out 
in  detail: 

fl  (xd)  =  m  (u  -  vr  +  wq) 

hj  =  Fx  -  2  Wj  sin  6 

T  .  ... 

Xd  =  ^U’  V>  w’  V’  r’ 

where  u,  v,  w,  u,  v,  w,  are  the  longitudinal,  lateral  and  normal  velocity  and  acceleration  components 
respectively; 

p,  q,  r,  p,  q,  r  are  the  roll,  pitch  and  yaw  angular  rates  and  angular  accelerations  respectively; 
the  second  term  in  hj  denotes  the  effect  of  the  weights  blown  from  the  ballast  tanks, 
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Fx  -  thrust  from  propeller  rpm 
m  -  mass  of  vehicle 
9  -  pitch  angle 

The  hydrodynamic  coefficient  vector  (cj)  for  this  equation  is  given  by  : 

T 

c>  "  (xqq’’xrr  xrp>  xu’  xvr  xwq>  xuu’  xw’  xww>  x5r5r>  x5s5s’  x5b5b^ 

where  5S,  5^  and  5r  denote  stern,  bow  and  rudder  surface  deflections.  The  Tj  vector  is  correspond¬ 
ingly  given  by 

r7  =  (q2 ,  r2 ,  rp,  u,  vr,  wq,  u2 ,  v2 ,  w2 ,  u2  5r2 ,  u2  5s2 ,  u2  5b2 ) 

The  remaining  equations  (i  =  2,  3,  4,  5,  6)  can  be  similarly  written  out.  It  is  not  necessary  explicitly 
to  represent  time  dependency  of  the  parameters  in  this  notation  since  time  can  always  be  introduced 
as  a  state  of  the  system. 
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V.  PROBLEM  FORMULATION 


The  utlimate  goal  is  to  determine  the  quality  and  quantity  of  instrumentation  required  to 
estimate  the  hydrodynamic  coefficient  vectors,  cj,  from  a  sequence  of  control  surface  deflections 
over  a  given  time  interval.  The  immediate  objective  of  this  work  is  to  develop  a  data-process- 
ing  method  that  will  aid  in  the  achievement  of  this  goal.  Note  that  the  control  inputs  (5)  are  initially 
specified  as  arbitrary.  It  is  of  course  desirable  that  they  be  sufficiently  random  such  that  the  signals 
which  form  the  coefficient  of  each  parameter  be  linearily  independent  of  the  signals  for  all  other  par¬ 
ameters  in  the  system.  An  additional  reason  for  randomness  of  the  control  sequence  will  be  explained 
later  in  relation  to  the  convergence  and  performance  of  the  proposed  coefficient  estimation  algorithm. 

The  system  states  which  define  the  motion  of  a  submarine  are  such  that  they  can  be  readily 
obtained  by  instrumentation;  e.g.,  gyros  and  accelerometers  can  be  placed  on  board  which  directly 
measure  u,  v,  w,  p,  q,  and  r.  Gravity  and  earth’s  spin  rate  are  of  course  subtracted  from  the  respec¬ 
tive  instrument  outputs.  If  the  initial  attitude  of  the  vehicle  is  established,  u,  v,  w  can  then  be 
obtained  by  what  is  equivalent  to  a  strapdown  inertial  navigation  mechanization,  while  p,  q  and  r 
can  be  directly  instrumented  or  obtained  from  the  gyro  outputs. 

The  essential  issue  here  is  that  estimates  of  the  dynamic  states  can  be  made  without 
imposing  (as  part  of  the  problem)  the  constraints  of  the  force-moment  equations.  Any  instrument 
employed  to  measure  a  variable  is  essentially  a  biased  estimator;  e.g.,  an  accelerometer  is  a  biased 
estimator  of  acceleration  or  a  gyro  is  a  biased  estimator  of  attitude  rate.  For  any  vehicle  employing 
inertial  navigation  the  standard  mechanizations  are  biased  estimators  of  the  dynamic  states. 
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VI.  ESTIMATION  ALGORITHM 


The  hydrodynamic  parameter  identification  problem  can  be  treated  by  the  methods  outlined 
in  Appendix  A.  The  technique  is  to  obtain  a  priori  estimates  of  the  dynamic  variables  (kinematic 
part  of  the  problem)  by  conventional  methods  (discussed  in  Section  X),  and  then  incorporate  the  con¬ 
straints  of  the  force-moment  equations  as  nonlinear  measurements  on  an  error  system  whose  states  consist 
of  vector  augmented  with  the  (Tj,  i  =  1,  2,  ...  6,  error  vectors.  The  resultant  error  state  vector 
x  is  now  given  by 


The  vector  x  will  denote  the  true  values  for  the  respective  elements. 


If  it  is  assumed  that  the  hydrodynamic  coefficients  are  constant  (cj  =  0),  then  the  error  dynamics 
is  reduced  to  specifying  the  function  for  x^.  The  simplest  case  is  the  assumption  of  noise-free, 


linear,  time-invariant  dynamics  for  xj, 
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This  is  valid  for  instances  where  the  duration  of  the  parameter  estimation  interval  is  short.  Where 
this  does  not  apply,  a  higher  degree  of  freedom  for  x^  is  obtained  by  making  u,  v,  w,  p,  q  and  7 
random  correlated  processes.  The  augmented  system  error  dynamics  can  now  be  summarized. 
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The  measurements  on  this  system  are  the  force-moment  equations  and  are  given  by, 


y.- 


ci  (xd,5)  +  hi-fi(xd),  i  =  1,2, 


where  yj  is  identically  zero.  When  estimated  quantities  for  q  and  xj  are  substituted  in  yj  the 
expressions  become  yp 


Therefore,  yj.  the  error  of  estimation  of  yj,  is  given  by, 


A/  A  ^ 

Yi  =  Vi-yi  =-yi- 

Hxpanding  y-  about  x;^  and  cj  and  setting  equal  to"yj, 

vi  =  ^rir$d,«),-jlxd-[fiUj).^d’JT  ~d  +  (?d,«)ci  (2) 

-*  +  higher  order  terms. 

where  denotes  vector  partial  differentiation  with  respect  to  x^.  These  expressions  can  be 
summarized  in  the  following  form: 


A 

y\ 


y2 


M  (x)  x  +  v 


(3) 


A 

y  6 


where  M(x)  is  a  6xn  matrix  which  is  a  function  of  x,  and  the  vector  v  is  introduced  as  white  measure¬ 
ment  noise.  Obviously  this  assumption  is  not  correct  —  but  is  necessary  in  order  to  apply  Kalman 
filtering  theory  .  Since  the  elements  of  M  are  functions  of  the  states  it  becomes  necessary  to  incor¬ 
porate  multiple  iterations  into  any  computational  algorithm.  A  major  practical  problem  is  the 
dimension  (n)  oFx!  The  number  of  coefficients  for  a  typical  vehicle  can  easily  exceed  100.  For  this 
reason  the  estimation  of  x  is  partitioned  into  six  separate  Kalman  filtering  problems  which  share  the 
covariance  matrix  for  x^, 


The  six  partitioned  error  state  systems  are  given  by, 


xd 

xd 

1 

1 

tsj 

Cj 

5 

~2 

)  '  *  * 

1 

lo 

_ i 

Let  Pj  denote  the  error  covariance  matrix  for  system  i.  A  computational  algorithm  which  successively 
applies  the  discrete  Kalman  filtering  equations  can  now  be  formulated.  It  will  be  assumed  that  the 
initial  estimates  for  x^  have  been  obtained  at  discrete  points  in  time  (t^)  and  are  stored  into  an  array. 

It  is  further  hypothesized  that  convergence  of  the  algorithm  may  be  sensitive  to  the  order  in  which  the 
error  systems  are  processed.  Therefore,  the  equation-measurement  processing  sequence  is  initially 
specified. 


In  developing  the  following  procedure  it  was  assumed  that  little  or  no  information  on  the  cj 
coefficients  is  available  other  than  their  initial  variance;  i.e.,  the  initial  estimates  for  cj  are  zero  and 
the  covariance  matrix  is  diagonal  and  its  elements  are  equal  to  the  initial  variances.  This  situation 
tends  to  invalidate  the  approximations  implicit  in  equation  3;  i.e.,  the  higher  order  terms  of 
equation  2  are  not  insignificant.  If  the  partial  derivative  terms  (coupling  the  x^  states)  are  initially 
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evaluated,  and  an  attempt  is  made  on  the  first  pass  accurately  to  balance  the  force  or  moment  mea¬ 
surement  equation,  then  divergence  of  the  extended  Kalman  filter  is  guaranteed.  To  prevent  this 
condition  from  occurring,  the  procedure  is  first  to  process  the  set  of  measurements  y(t]<)  with  a 
large  value  for  Qj,  where  Qj  =E  j  v{  U[  }  ,  and  on  successive  passes  over  the  same  equations  introduce 
the  state  terms  by  evaluating  the  respective  partial  derivatives.  On  these  additional  passes  the 
value  of  Qj  is  reduced.  Qj  is  the  noise  covariance  matrix,  and  in  this  case  is  a  scaler. 

The  important  issue  is  to  obtain  from  Pj  a  matrix  which  accurately  reflects  the  state 

of  uncertainty  in  The  resulting  estimates  for  x^  obtained  from  a  particular  equation  are  then 
used  to  update  the  values  for  x^.  The  starting  covariance  for  ^x^  in  the  next  system  Pj,  i  4=-  j,  is 
obtained  from  Pj.  For  the  deterministic  error  dynamics  postulated  by  equation  1 ,  transferring  the 
terminal  P^  to  initial  time  presents  no  problem;  but  for  a  more  general  stochastic  error  structure 
it  is  necessary  to  construct  multiple  point  smoothing  into  the  algorithm.  The  remaining  equation- 
measurements  are  processed  in  a  similar  manner.  The  primary  motivation  for  this  procedure  is 
that  the  six  equations  of  motion,  processed  sequentially,  can  provide  sufficient  observability  for  x^ 
so  that  improved  estimates  of  cj  are  obtained  if  the  force-moment  equations  are  reprocessed. 

The  basic  algorithm  then  is  a  sequence  of  multiple  sweeps  on  six  error  state  systems  where 
the  measurement  on  each  system  consist  of  the  corresponding  force  or  moment  equation.  After 
each  of  the  six  systems  has  been  multiple-iterated  once,  the  computations  are  repeated  with  fewer 
individual  sweeps  since  the  greatest  convergence  will  occur  during  the  initial  sweeps.  From  simula¬ 
tions  of  this  algorithm  it  was  established  that  only  one  additional  sweep  through  each  system  was 
required  to  obtain  the  parameter  vectors  Cj. 

Convergence  of  this  technique  is  directly  related  to  obtaining  unbiased  estimates  for  Cj  on 
the  first  pass,  where  the  value  of  Qj  is  large.  From  this  it  is  asserted  that  the  sequence  of  control 
surface  deflections  should  be  as  random  as  possible  in  an  attempt  to  whiten  the  higher  order  terms 
of  equation  2,  thereby  making  valid  the  assumptions  of  the  discrete  Kalman  filter. 

The  major  source  of  instability  in  this  algorithm  is  the  process  of  updating  estimates  with 
the  successive  sweeps.  This  may  account  for  the  instability  noted  at  present  in  all  applications  of 
extended  Kalman  filtering.  It  was  discovered  that  false-observability  was  created  if  the  corrections 
to  the  Cj  obtained  from  the  current  sweep  were  used  in  that  sweep  to  compute  the  partials  coupling 
the  state  into  the  measurement.  This  is  obvious  when  it  is  considered  that  the  Cj  vectors  are 
constant;  to  introduce  the  convergence  dynamics  of  the  Cj  vectors  into  the  measurement  structure 
is  a  gross  modeling  error.  This  may  be  the  reason  why  some  researchers  have  concluded  that  the 
extended  Kalman  filter  is  unstable  or  very  sensitive;  and  have  resorted  to  least-squares  methods 
which  inherently  avoid  this  observability  distortion. 

This  difficulty  is  avoided  by  simply  using  the  end  results  from  the  previous  sweep  in 
evaluating  the  partials  on  the  current  sweep. 

The  logical  flow  diagram  presented  as  Figure  1  shows  the  interrelationships  of  the  various 
sweeps,  iterations  and  updating  procedures. 

A  separate  subroutine  (YWIG)  was  constructed  into  the  program  for  the  computation  of  yj 
(Figure  1).  This  made  the  computation  of  the  partials  (numerically)  convenient  by  calling  YWIG  for 
two  different  values  of  a  particular  dynamic  variable  -  differencing  the  result  and  dividing  this  by  the 
difference  between  the  two  values  of  the  dynamic  variable. 


1 1 
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Figure  1.  Flow  Diagram  for  Hydrodynamic  Coefficient  Estimation  Algorithm 

12 


VII.  INTRODUCTION  TO  SIMULATIONS 


The  approach  taken  by  BAC  to  evaluate  this  algorithm  was  to  initially  solve  the  simplest  pro¬ 
blem  possible.  This  is  the  case  where  all  the  states  x^  are  constant,  i.e.,  F=0.  The  main  reason  for 
doing  this  is  that  the  time  varying  problem  would  require  greater  programming  and  larger  execution 
time.  It  was  felt  that  under  the  limited  funding,  a  realistic  and  useful  goal  would  first  be  established. 
Also,  if  a  simplified  version  of  the  problem  could  not  be  solved  then  it  would  be  economically  waste¬ 
ful  to  expend  effort  developing  sophisticated  processing  algorithms  which  had  no  chance  of  success. 
This  simplified  problem  may  turn  out  to  be  the  best  way  in  which  the  real  world  submarine  identifica¬ 
tion  problem  can  be  approached.  The  real  situation  is  often  characterized  by  the  occurrence  of  events 
that  are  unexpected,  i.e.,  no  matter  what  the  assumed  mathematical  model  for  the  submarine  is,  there 
is  a  significant  probability  that  it  is  not  quite  correct.  If  some  hydrodynamic  coefficients  are  not 
constant,  how  then  should  they  be  modeled?  Are  they  functions  of  depth,  speed,  etc.?  Higher 
degrees  of  freedom  can  be  built  into  an  algorithm  by  assuming  random  correlated  (stochastic)  model¬ 
ing  for  any  coefficient(s);  this  places  a  much  greater  decision  making  requirement  on  the  program 
operator  as  to  where  a  higher  degree  of  freedom  is  needed.  This  is  a  difficult  area  because  it  is  con¬ 
ceivable  that  at  some  level  of  freedom  nothing  is  gained,  i.e.,  submarine  coefficients  can  be  obtained 
to  fit  any  set  of  random  signals. 

These  considerations  make  the  method  outlined  above  attractive  as  a  means  of  attacking  the 
real  world  problem.  It  follows  that  it  would  be  more  useful  and  practical  to  break  up  an  hour  of 
submarine  trajectory  data  into  short  separate  intervals  and  process  the  data  as  though  it  was  from 
different  submarines,  rather  than  have  a  complex  time-varying  identification  algorithm  to  process 
the  data  as  taken  from  one  submarine.  This  general  approach  had  been  taken  in  reducing  flight 
data  on  BAC  Hipernas  inertial  navigators  and  was  found  to  be  the  most  beneficial.  To  this  goal  then, 
it  is  hoped  that  the  following  simulations  will  contribute  significantly. 
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VIII.  SIMULATIONS 


The  following  are  results  from  simulations  of  a  typical  submerged  vehicle  with  zero  initial 
conditions  —  except  for  u  =  50  ft/sec.  The  control  surface  commands  are  basic  pseudo  random 
sequences  with  a  5  sec  pulse  interval.  The  thrust  is  such  that  u  is  decreased  to  20  ft/sec  over  a 
10-minute  flight  time.  Salient  features  of  the  resultant  trajectory  are  shown  in  Figure  2.  The 
sequential  Kalman  filter  computation  interval  is  six  seconds;  processing  of  data  was  performed  for 
a  five  minute  flight.  It  required  five  minutes  of  computing  time  on  an  IBM  360/65  to  compute  a 
trajectory  and  estimate  the  parameters  and  bias  values.  The  number  of  states  for  each  error  system 
and  the  order  in  which  they  were  processed  are: 


Axial  Force  System 

24 

Rolling  Moment  System 

25 

Pitching  Moment  System 

28 

Normal  Force  System 

28 

Lateral  Force  System 

27 

Yaw  Moment  System 

26 

These  values  include  the  1 2  dynamic  states,  x^.  It  might  be  argued  that  the  order  of  these  separate 
systems  should  be  lower  since,  for  each  system,  some  of  the  x^  variables  are  not  directly  involved. 
However,  the  algorithm  advanced  here  does  couple  the  errors  of  these  variables  (in  the  error  covariance 
matrix)  as  the  systems  are  sequentially  processed  and  the  covariance,  for  the  x  from  processing  a 
given  system  is  used  to  start  the  processing  of  the  following  system. 

The  dynamic  states  (x^)  were  generated  for  the  true  coefficients  with  the  submarine  simula¬ 
tion  program  and  submarine  parameter  data  supplied  to  BAC  by  NSRDC  [3] .  The  hydrodynamic 
coefficient  estimation  program  was  written  to  operate  on  a  two  dimension  array  that  consist  of 
sampled  values  of  the  following  dynamic  variables  at  discrete  points  in  time, 

u,  v,  w,  u,  v,  w,  p,  q,  r,  p,  q,  r,  5r,  6s,  5  b, 

PRMNI,  FX,  FY,  FZ,  FM,  FN,  FK 

The  first  1 2  variables  are  the  dynamic  variables  (xj)  that  define  the  motion  of  the  vehicle;  6 r,  5  s,  5  b 
are  the  control  surface  deflections;  PRMNI  is  (n'  -  1 )  where  n'  is  the  ratio  of  ordered  speed  to  current 
axial  velocity.  In  the  set  of  hydrodynamic  coefficients  supplied  to  BAC,  the  parameters  associated 
with  PRMNI  were  zero,  therefore  at  this  time  estimation  of  these  terms  have  not  been  investigated. 

The  remaining  terms  in  the  above  array  are  defined  as  follows: 


FX  = 

-  2W-  sin  9  +  FXP 

FY  = 

2Wj  sin  0  cos  9  +  FYS 

FZ  = 

ZWj  cos  9  cos  0  +  FZS 

FM  = 

B(Zg)  sin  9  -  SWjXjj  cos  9  cos  0  +  QYS 

FN  = 

iWjXyj  cos  9  sin  0  +  QZS 

FK  = 

B(Zg)  cos  9  sin  0  +  QXS 
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Time  in  Seconds 

Figure  2.  Trajectory  Data  and  Control  Surface  Deflections  for 
Parameter  Estimation  Simulations 


The  S  terms  here  are  sea  state  effects  which  the  simulation  program  sets  equal  to  zero;  9  and 
0  are  pitch  and  roll  of  the  submarine.  As  a  first  step  in  the  development  of  the  program  it  was  felt 
that  identification  of  the  parameters  should  be  investigated  with  knowledge  of  Zg,  0  ,  <j>,  FXP  and 
2Wj  (sum  of  weights  blown).  It  is  recognized  that  any  real  world  situation  has  errors  in  these  terms  - 
also  errors  in  knowledge  of  the  CG  exist.  The  program  is  structured  such  that  extending  the  degrees 
of  freedom  to  include  these  errors  is  a  trivial  matter. 

Before  the  data  outlined  above  is  reduced  (estimating  the  hydrodynamic  coefficients)  the 
variables  are  first  corrupted  by  bias  and  random  errors. 

The  results  presented  in  tables  1  through  7  are  for  one  set  of  bias  errors;  these  are  given  in 
Table  7.  Table  7  also  shows  how  well  the  bias  errors  were  estimated  (by  the  algorithm)  for  two 
levels  of  white  noise,  (1)  zero  and  (2)  as  follows: 


Au  = 

0.01  ft/sec 

A5s  = 

0.0005  rad 

AV  = 

0.005  ft/sec 

Au  = 

0.0003  ft/sec2 

Aw  = 

0.005  ft/sec 

Av  = 

0.0003  ft/sec2 

Ap  = 

0.5(1  OF5  rad/sec 

Aw  = 

0.0003  ft/sec2 

Aq  = 

0.5(1 0FS  rad/sec 

Ap  = 

0.2(1 0FS  rad/sec 

Ar  = 

0.5(1 0FS  rad/sec 

Aq  = 

0.2(1 0FS  rad/sec 

A6r  = 

0.0005  rad 

Ar  = 

0.2(1 0FS  rad/sec 

A5b  = 

0.0005  rad 

These  random  errors  were  added  to  the  trajectory  variables  as  independent-uniformly  distributed 
processes.  Tables  1  through  6  summarize  the  accuracies  obtained  for  the  dimensionless  hydrodynamic 
coefficients  associated  with  each  force/moment  equation.  Those  hydrodynamic  parameters  that  do 
not  appear  in  the  table  were  given  to  BAC  by  NSRDC  as  zero  and  therefore  not  made  a  part  of  this 
study.  The  program  is  coded  for  the  n  '  coefficients  and  additional  ones  can  readily  be  added. 

The  convergence  dynamics  for  every  coefficient  estimated  is  of  interest  because  it  readily 
identifies  the  observability  structure  of  the  parameters.  This  information  can  be  used  to  ascertain 
whether  or  not  a  particular  coefficient  would  be  observable  if  it  were  not  constant  over  the  estimation 
time  interval.  For  this  reason  normalized  plots  of  the  convergence  to  each  coefficient  is  provided  at 
the  end  of  this  report.  Of  particular  significance  is  the  convergence  sequence  depicted  in  the  plots 
and  the  corresponding  noise  levels,  Qj,  which  control  the  gains  in  the  Kalman  filter.  The  form  of  Qj 
and  the  operational  details  of  the  computer  algorithms  looping  pattern  are  summarized  below. 

1 .  Each  system-equation  was  first  iterated  five  times  in  the  following  equation  number 
sequence,  1-4-5-3-2-6.  The  end  covariance  from  the  first  sweep  was  stored  and  used  for  initializing 
the  succeeding  sweeps. 

2.  Estimates  of  x^,  which  are  obtained  starting  with  the  second  sweep,  are  reset  to  zero 
before  starting  the  succeeding  sweeps.  After  the  final  sweep  on  an  equation-system,  the  estimates 
are  used  to  update  the  array  of  xcj(tj<).  This  array  is,  at  the  start  of  the  algorithm,  the  sum  of  the 
true  dynamic  variables  (x^)  plus  the  bias  and  random  elements  described  above. 

3.  As  in  estimating  the  x^  variables,  corrections  to  the  hydrodynamic  coefficients  are 
obtained  from  any  sweep,  and  at  the  final  time  of  each  sweep  the  estimates  are  used  to  update  the 
respective  coefficient  array.  This  is  necessary  to  avoid  introducing  the  convergence  dynamics  into 
the  measurement  structure. 
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TABLE  1 

SUMMARY  OF  AXIAL  FORCE 
DIMENSIONLESS  COEFFICIENT  ESTIMATION 


Estimated  Values 

With  Bias 

%  Estimated 

(Based  Upon  Covariance) 

Parameter 

True  Value 

No  White  Noise 

White  Noise 

No  White  Noise 

With  White  Noise 

xqq 

-1.7(10)"4 

-1.695(1  Or4 

-1 .6 1 9(  1 0)-4 

96.23 

56.23 

xrr 

-8.6(10f5 

-9.553(1 0)"s 

-2.737(10F5 

47.42* 

1.45 

Xrp 

2.7(10f4 

2.702(10)'4 

2.744(1  Of4 

98.71 

83.09 

-1.42(1  or4 

-1.392(1  or4 

0.1 562(1 0)~4 

91.94 

29.94 

Xvr 

i.id  or2 

1.093(1  or2 

1.072(1 0)"2 

97.38* 

91.56* 

XWq 

-7.46(1  or3 

-7.455(1  Of3 

-7.209(1  Of3 

99.70 

97.05 

Xuu 

0. 

- 

- 

- 

- 

Xw 

6.5(1  or3 

6.388(1 0)"3 

6.025(1 0)-3 

92.58* 

60.42* 

Xww 

1.92(1  or3 

1.929(1  or3 

2. 174(1  or3 

99.02 

89.56 

1 

•2.89(10)"3 

-2.887(1 0)“3 

-2.826(1 0)-3 

99.53 

94.28 

-2.46(1  or3 

-2.458(1 0)"3 

.  -2.428(1 0)"3 

99.27 

90.34 

X§b5b 

-2.59(1  or3 

-2.585(1 0)“3 

-2.370(1  Of3 

99.49 

93.43 

TABLE  2 

SUMMARY  OF  LATERAL  FORCE 
DIMENSIONLESS  COEFFICIENT  ESTIMATION 


Parameter 


True  Value 


•1.1(1  Of2 
-7.3(1  or3 
7.46(1  or3 
3.57(1  or3 
-1.30  or3 
-6.74(1  or4 
0. 

-0.06 

-0.0213 

-0.00652 

-0.066 

0. 

0. 

-0.02075 


Estimated  Values 

With  Bias 

%  Estimated) 

(Based  Upon  Covariance) 

No  White  Noise 

White  Noise 

No  White  Noise  With  White  Noise 

1.694(1  or4 

0.298(1  or4 

93.68* 

56.65 

-2.700(1  or4 

-2. 584(1  OF4 

99.82 

97.86 

8.539(1 0)"s 

-7 .903(1 0)_s 

87.85* 

24.61 

•I.IOO(IO)"2 

-0.995(1  or2 

99.67 

96.06 

-7.299(1 0)"3 

-6 .464(1  OF3 

99.32 

91.67 

7.457(1 0)"3 

6.6 13(1  or3 

99.61 

95.29 

3.569(1 0)“3 

3.227(  10)-3 

98.93* 

92.48 

-1.299(1 0)"3 

-0.9307(1  or3 

98.69 

88.22 

-6.740(  10)"4 

-6.832(1  Of4 

99.62 

96.86 

-0.05998 

-0.05271 

99.46 

95.64 

-0.02130 

-0.02161 

99.56 

96.34 

-0.006517 

-0.005578 

99.38 

95.84 

-0.06598 

-0.05919 

99.66 

96.81 

— 

— 

—  — 
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TABLE  3 

SUMMARY  OF  NORMAL  FORCE 
DIMENSIONLESS  COEFFICIENT  ESTIMATION 


Estimated  Values 

With  Bias 

%  Estimated 

(Based  Upon  Covariance) 

Parameter 

True  Value 

No  White  Noise 

White  Noise 

No  White  Noise 

With  White  Noise 

Zrr 

-0.00178 

-0.001782 

-0.001794 

96.24* 

82.47* 

Zrp 

-5.33(10)-* 

-5.332(  10)-4 

-3. 164(1  or4 

97.98 

77.83 

Zq 

•1.7(1  or4 

-1.699(1  or4 

-1.559(10)  "4 

98.99 

88.94 

Zw 

•0.00746 

-0.00746 

-0.007418 

99.96 

99.56 

zvr 

0.00203 

0.002015 

0.00151 

78.39* 

7.71* 

2wjq| 

0.0062 

0.0062 

0.006356 

99.69 

96.22 

Zq 

-0.0044 

-0.0044 

-0.004364 

99.87 

98.65 

Z|q|6s 

0. 

- 

- 

- 

- 

Zvp 

-0.0068 

-0.0068 

-0.006043 

99.46 

94.62 

zAvp 

0. 

- 

- 

- 

- 

N 

c 

c 

-1.1 1(1 0)-4 

-1.11(1  or4 

-1.121(1  or4 . 

99.85 

98.00 

N 

< 

< 

0.066 

0.06598 

0.06442 

98.86* 

94.68 

Zw|w| 

-0.03184 

-0.03184 

-0.03165 

99.83 

97.95 

zw 

-0.01059 

-0.01059 

-0.0106 

99.93 

99.24 

Zfis 

-0.00542 

-0.00542 

-0.00534 

99.91 

98.95 

Z5b 

-0.00267 

-0.00267 

-0.002665 

99.95 

99.39 

Z|w| 

0. 

- 

- 

- 

- 

ZWw 

0. 

- 

— 

— 

— 

zvs 

-0.02075 
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TABLE  4. 

SUMMARY  OF  ROLL  MOMENT 
DIMENSIONLESS  COEFFICIENT  ESTIMATION 


Estimated  Values 

With  Bias 

% 

Estimation 

(Based  Upon  Covariance) 

Parameter 

True 

Value 

No  White 

Noise 

White 

Noise 

No  White 

Noise 

With  White 

Noise 

kp!pI 

-3.4(1  or6 

-3.398(1  Of6 

-3.828(1 0)'6 

98.72 

78.14 

Kqr 

-1.020  or4 

-1.02(1  or4 

-1.028(1  Of4 

99.93 

98.90 

Kp 

-3.3(1  or6 

-3.3(1  or6 

-3.194(10)'6 

99.90 

98.29 

Kf 

-1.44(1  or5 

-1.44(1  or5 

-1.487(1  or5 

99.88 

97.94 

Kv 

-2.7(1  or4 

-2.70(1  or4 

-2.623(1  Of4 

99.91 

98.49 

Kwp 

2.7(10)"4 

2.7(10)-4 

2.629(1 0)"4 

99.92 

98.66 

Kr 

0. 

- 

- 

- 

- 

Kp 

-3.6(1  or5 

-3.6(1  or5 

-3. 582(1 0)'s 

99.96 

99.41 

Kuu 

0. 

- 

- 

- 

- 

Kv|v| 

-1.037(1  or3 

-1.037(1  or3 

-0.9744(1  Of3 

99.85 

97.50 

Kv 

-3.87(1  Of4 

-3.871(1  or4 

-3.845(1 0)'4 

99.96 

99.28 

Kfir 

6.4(1  or5 

6.4(1  or5 

5.993(1 0)"s 

99.79 

96.45 

Kwv 

3.54(1  or4 

3.539(1  Of4 

3.931  (10)"4 

99.63 

93.72 
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TABLE  5 

SUMMARY  OF  PITCH  MOMENT  DIMENSIONLESS 
COEFFICIENT  ESTIMATION 


Parameter 


Mrr 

Mrn 


Mq 

MlqISs 

Mwn 


-0.00112 
9.92(1 0)' 5 
-3.9(1  Of4 
0 

-1.73(10)"4 

-0.00362 

-0.0022 

-0.00255 

0 

8.44(1 0)'4 
4.2(1  or5 
0.0146 
-0.00473 
0.0032 
-0.00258 
5.2(1  or4 
0 
0 


%  Estimation 

Estimated  Values 

(Based  Upon 

With  Bias 

Covariance) 

No  White 

White 

No  White 

With  White 

Noise 

Noise 

Noise 

Noise 

-0.001121 

-0.0009485 

99.92 

77.70* 

9.944(1  Of5 

13.1  k  iors 

99.87 

58.62 

-3.9(1  or4 

-s^nor4 

99.99 

98.17 

-1.73(10f4 

-1.645(1 or4 

99.97 

94.18 

-0.003626 

-0.003081 

99.84 

55.19* 

-0.0022 

-0.002241 

99.98 

94.79 

-0.00255 

-0.002504 

99.99 

99.16 

8.447(1  Of4 

10.50(1  or4 

99.94 

84.32 

4.2(1 0)_s 

4.275(1 0)-s 

99.96 

98.13 

0.01459 

0.01427 

99.93 

80.79* 

-0.004732 

0.005042 

99.98 

94.88 

0.0032 

0.003177 

99.99 

99.10 

-0.00258 

-0.002529 

99.99 

99.13 

5.2(1  or4 

5.086(1  Of4 

99.99 

98.97 

- 

- 

- 

- 

TABLE  6 

SUMMARY  OF  YAW  MOMENT 
DIMENSIONLESS  COEFFICIENT  ESTIMATION 


Estimated  Values 

With  Bias 

%  Estimated 
(Based  Upon 

Covariance) 

Parameter 

True 

Value 

No  White 

Noise 

White 

Noise 

No  White 

Noise 

With  White 

Noise 

Npq 

-^(ibr4 

-3.872(10)"4 

-o.sssoor4 

99.95 

89.92* 

Np 

-7.2(1  or6 

-7. 184(1  Of6 

-6.59(1  Of6 

99.92 

83.09 

Nr 

-4 .92(1 0)-4 

-4.922  (10F4 

-1.426(10^ 

99.95 

91.60* 

N; 

3.12(10)'4 

3.131  (10F4 

0.888(1  OF4 

99.86 

71.55 

NlJr 

-4 .45(1  OF3 

-4.450(1  Of3 

-3,073(1  Of3 

99.97 

93.42 

Nwp 

-1 .73(1  or4 

-1.739(10]^ 

-0.424  (10)"4 

99.79 

57.13 

Nr 

-3.2(1  or3 

-3.200(10f3 

-1.986(1  Of3 

99.97 

95.50* 

N|rl5r 

1.3(1  or3 

1.300(10)“3 

0.81 9(1 0)-3 

99.97 

95.01 

Np 

-2.2(1or4 

-2.200(10)~4 

-1.445(10)~4 

99.98 

95.50* 

Nuu 

0. 

- 

- 

- 

0.0135 

0.01351 

0.00681 

99.96 

92.66 

Nv 

-0.00713 

-0.007130 

-0.00448 

99.97 

95.30* 

Nfir 

-0.00333 

-0.003331 

-0.00204 

99.97 

95.45* 

Nwv 

0.0146 

0.01461 

0.00901 

99.97 

95.01 

TABLE  7 

ESTIMATION  OF  BIAS  ERRORS 


Estimated  Values 

With 

%  Estimated 
(Based  Upon 

Covariance) 

True 

No  White 

White 

No  White 

With  White 

Parameter 

Bias  Value 

Noise 

Noise 

Noise 

Noise 

ft/sec 

3.0 

3.000 

2.939 

99.98 

Av  ft/sec 

0.3 

0.3001 

0.303 

99.99 

Aw  ft/sec 

0.3 

0.3002 

0.2994 

99.96 

96.14 

Au  ft/sec2 

0.0032 

0.00316 

0.00203 

97.48 

70.4 

Ay  ft/sec2 

0.0032 

0.00320 

0.00318 

99.58 

92.43 

Aw  ft /sec  2 

0.0032 

0.00318 

0.00287 

96.84 

57.50 

Ap  rad/sec 

0.5(1  Of4 

0.504(10f4 

0.237(10)~4 

98.39 

36.68 

Aq  rad/sec 

0.5(10r4 

o^nor4 

-0.066(1  or4 

95.89 

20.47 

Ar  rad /sec 

0.5(1 0)"4 

0.496  (10)"4 

0.295(1  OF4 

99.65 

82.4 

Ap  rad/sec2 

0.2(1  or5 

0.068(10r5 

-0.007(1 0)_S 

45.20 

0.11 

Aq  rad/sec2 

o.2(iors 

0.237(1  Of5 

-0.006(1  or6 

79.26 

0.40 

Af  rad /sec  2 

o.2(ior5 

0.205(10f5 

-0.055(10rS 

97.03 

3.06 
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4.  Initialization  of  the  covariance  matrix  for  a  given  system-equation  consist  of  diagonaliz¬ 
ing  the  covariance  matrix  for  q  and  zeroing  the  off-diagonal  elements  between  x^  and  q.  This  is 
done  only  for  the  first  sweep;  succeeding  sweeps  are  initialized  by  setting  the  covariance  of  the 
equation-system  equal  to  that  obtained  from  the  end  of  the  initial  sweep  (where  Qj  was  high).  The 
covariance  for  x^  is  initialized  as  diagonal  and  is  not  altered  when  sweeps  on  a  new  equation-system 
are  initiated. 

5.  The  measurement  noise  levels  on  the  initial  five  sweeps  were  computed  as  follows: 

Q-  =  R^2  (SR-)2  for  first  sweep 

Qj  =  /Rj  +  CRjSRjHl  ’  fori  =  2,3,4 

where  I  is  the  sweep  number  and  SRj  is  a  dimensionless  gain  factor  introduced  to  set  the  level  of  Qj. 
The  fifth  sweep  is  made  adaptive  with  respect  to  noise  by  computing  Qj  from  an  inputed  array 
reflecting  the  expected  values  of  the  white  noise  errors  in  x^,  and  the  sensitivity  functions  coupling 
the  bias  errors  into  the  measurement.  The  selection  of  Qj  sequence  is  somewhat  aribtrary;  the  partic¬ 
ular  choice  made  here  is  based  upon  the  intuitive  notion  that  the  algorithm  should  converge  (most 
naturally)  on  a  parabolic  type  trajectory. 

6.  After  processing  each  system  five  times,  one  additional  sweep  is  made  through  each 
system  with  Qj  computed  as  in  the  final  loop  of  the  first  sequence  of  sweeps.  The  values  of  Rj  and 
SRj  used  in  the  simulations  are  given  below, 

No  White  Noise  White  Noise 

R!  =  3001b  ,  SRi  =  50  Rj  =  3001b  ,  SRX  =  50 

R2  =  300  1b  ,  SR2  =  50  R2  =1000  lb  ,  SR2  =  200 

R3  =  300  lb  ,  SR3  =50  R3  =  600  lb  ,  SR3  =  200 

R4  =  3000  ft-lb,  SR4  =  50  R4  =30,000  lb,  SR4  =  50 

Rs  =  3000  ft-lb,  SRs  =100  Rs  =30,000  lb,  SRs  =  100 

R6  =  3000  ft-lb,  SR6  =100  R6  =100,000  ft-lb,  SR6  =  100 

It  should  be  mentioned  that  these  values  were  arrived  at  by  a  closed  loop  operation 
between  the  computer  program  developer  and  the  computer  output.  The  degree  to  which  this  is 
necessary  is  directly  related  to  the  degree  of  ignorance  in  the  apriori  estimates  of  the  hydrodynamic 
coefficients.  Instability  or  a  biased  performance  of  the  algorithms  will  be  experienced  if  these 
values  are  underestimated.  This  occurs  for  two  reasons:  (1 )  the  white  noise  corrupted  dynamic 
variables  are  used  in  constructing  the  measurement  matrix,  and  if  Qj  is  too  low  false  observability 
will  be  generated  in  the  covariance  matrix,  (2)  the  measurement  error,  y j,  cannot  be  forced  to  zero 
or  to  the  level  of  white  noise  without  being  able  to  compute  reasonably  accurate  partial  derivatives 
that  couple  the  bias  errors  into  the  equation-measurement.  It  should  be  emphasized  that  considerable 
latitude  does  exist  in  the  selection  Rj  and  SRj,  and  as  a  general  rule  it  is  safe  to  select  them  on  the 
high  side.  This  flexibility  results  from  the  form  of  Qj,  which  varies  Qj  through  a  large  dynamic 
range,  and  from  the  adaptive  feature  described  in  step  5  above  that  optimizes  the  value  of  Qj  on  the 
fifth  sweep,  and  when  the  equation-systems  are  reprocessed  on  the  final  pass. 

It  was  mentioned  in  step  4  above  that  the  starting  covariance  matrix  for  q  is  diagonal; 
for  this  study  the  value  of  these  diagonal  elements  were  obtained  by  squaring  twice  the  value  of  the 
hydrodynamic  coefficients  supplied  to  BAC  by  NSRDC.  Correspondingly,  the  initial  estimates  of  all 
coefficients  were  set  at  zero.  This  was  done  to  establish  a  very  conservative  apriori  knowledge  require¬ 
ment  for  starting  the  algorithm. 
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IX.  DISCUSSION  OF  RESULTS 


The  results  in  Tables  1-7  for  the  no  white  noise  case  were  obtained  with  the  gains  of  the 
Kalman  filter  (on  the  final  sweep)  computed  for  a  measurement  noise  level  of  100  lb  for  the  force 
equations  and  500  ft-lb  for  the  moment  equation.  Therefore,  the  final  covariance  matrix  reflects 
this  level  of  error  and  the  estimates  will  be  better  than  the  columns  designated  “%  Estimated”  (in 
the  tables)  indicate.  These  columns  are  computed  as  follows: 


%  Estimated  = 


/,  -EH 

y  V  pjj  ( lo ) 


(100) 


where  Pj  j  (  )  denotes  the  respective  diagonal  element  of  the  covariance  matrix  for  a  coefficient 

at  the  time  indicated. 


The  degree  to  which  a  particular  coefficient  can  be  estimated  (observability)  is  related  to  two 
physical  situations;  the  effect  of  a  coefficient  in  the  equation-measurement  may  be  (1  )  equivalent  to 
a  linear  combination  of  other  coefficients  or  (2)  at  a  level  that  is  within  the  measurement  noise.  It 
is  of  interest  then  to  determine  the  nature  of  an  observability  problem  associated  with  any  parameter. 
This  information  is  available  from  the  off-diagonal  elements  of  the  final  covariance  matrix.  The 
correlation  between  two  variables  are  obtained  from  the  following  expression, 


cor 


ij 


^  p«  pjj 


where  Pjj  is  the  i'  th  row  and  j  '  th  column  off-diagonal  element  of  the  covariance  matrix.  Variables 
that  are  in  combination  with  a  correlation  greater  than  0.9  are  designated  by  an  asterisk  (  *  )  in 
Tables  1  through  6.  For  coefficients  where  the  “%  Estimated”  does  not  indicated  good  estimation, 
and  the  variable  is  not  correlated  with  another  variable,  it  can  be  inferred  that  this  coefficient  does 
not  contribute  significantly  to  the  submarine  trajectory. 


A  situation  that  can  occur  is  that  some  variables  are  redundant  unless  the  level  of  random 
error  is  zero  or  very  small.  This  is  illustrated  by  the  pitch  equation  results  for  parameters  Mrr,  Mvr 
and  Mvy.  A  similar  effect  between  variables  Xrr,  Xvr  and  Xw  exist  for  the  axial  force  equations. 
But  in  this  case  the  correlation  for  Xrr  did  not  develop  in  the  white  noise  case.  This  indicates  that 
the  influence  of  the  Xrr  term  on  the  submarine  is  very  small.  To  obtain  a  quantitative  measure  of 
the  influence  of  these  terms,  Table  8  gives  the  nominal  influence  of  the  white  noise  variables  on  the 
balancing  of  each  equation.  The  term  nominal  is  used  here  because  the  effect  is  of  course  variable 
along  the  submarine  trajectory. 

In  Tables  2  and  3,  Yvs  and  Zvs  are  sinusoidal  hydrodynamic  effects  at  the  vortex  shedding 
frequence  (cc).  The  value  of  co  supplied  to  BAC  by  NSRDC  was  zero,  therefore  these  term  did  not 
excite  the  system  and  consequently  are  not  estimated. 


*  In  effect,  an  asterisk  indicates  that  the  functions  which  form  the  coefficients  of  the  marked 
parameters  become  correlated  at  the  level  of  the  parameter  value  indicated,  thus  preventing 
further  convergence. 
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TABLE  8 

EFFECT  OF  WHITE  NOISE  ERRORS 
ON  DYNAMIC  EQUATIONS 


Equation 

Effect 

Effect  of  Errors 

Number 

of  Xj  (White) 

in  Control  (White) 

1  (Axial) 

190  lbs 

1200  lbs 

2  (Lateral) 

900  lbs 

900  lbs 

3  (Normal) 

700  lbs 

900  lbs 

4  (Roll) 

750  lbs 

3500  ft  -  lbs 

5  (Pitch) 

50,000  ft  •  lbs 

150,000  ft -lbs 

6  (Yaw) 

80,000  ft  -  lbs 

150,000  ft-  lbs 

An  important  issue  of  submarine  parameter  identification  is  featured  in  the  results  (white 
noise  case)  for  the  yaw  coefficients.  It  is  noted  that  some  coefficients  are  not  estimated  within  the 
expected  accuracy  indicated  (Table  6).  The  question  immediately  arises  is  the  algorithm  failing  for 
the  level  of  noise  employed,  or  is  the  submarine  yaw  equation  non-unique.  The  test  was  to  simulate 
the  submarine  with  the  incorrect  coefficients  and  compare  the  resulting  trajectory  variables  with  the 
true  ones.  This  was  done  and  the  two  trajectories  matched  to  within  0.1  percent;  thus  indicating  that 
the  submarine  was  non-unique.  It  is  ascertained  that  the  reason  for  this  occurrence  is  related  to  the 
requirement  the  submarine  must  blow  weight  to  allow  a  unique  solution  for  the  yaw  equation.  This 
requirement  is  also  true  for  the  lateral  and  normal  force  equations.  The  basic  problem  is  that  the  in¬ 
ertia  terms  are  similar  to  some  hydrodynamic  coefficients,  e.g.,  the  yaw  equation  is 


I7  r  +  (I 


y  -  Ix)  PQ  ■  Swj  Xt.  cos  6  sin  <p 

j  j 

=  (-  2  5  N;)  r  + 

2  r  2 


N  pq  )  P<1 


+  other  coefficients 


Without  blowing  weight,  estimation  of  the  coefficients  must  come  from  the  information  con¬ 
tained  in  the  left  hand  side  of  this  equation.  But  the  left  hand  side  is  embedded  into  the  coefficient 
side  since  r  and  pq  appear  on  both  sides  thereby  unsealing  the  equation. 

Therefore  the  submarine  must  blow  weight  in  order  to  be  unique.  What  is  happening,  in  the 
results  of  Table  6,  is  that  the  selected  level  of  noise  is  starting  to  mask  the  equation’s  observability 
for  the  weight-blowing  profile  depicted  in  Figure  1.  This  situation  in  the  equations  also  forced  the 
computer  program  to  be  executed  in  double  precision.  Note  also  that  the  processing  sequence  for  the 
equation-systems  list  first  those  equations  that  do  not  require  mass  changes  to  be  observable.  In¬ 
stability  of  the  algorithm  was  experienced  for  sequences  that  processed  the  lateral  and  normal  force 
equation  before  the  pitch  and  roll  equations. 

This  relationship  between  uniqueness,  weight  blown  and  noise  level  is  an  area  requiring  further 
investigation  along  with  those  listed  in  the  introduction  to  simulations  section.  Another  simulation 
was  made  with  the  white  noise  a  factor  of  ten  higher  than  those  listed  previously.  In  this  case  estimates 
for  Equation  2  were  degraded  along  with  Equation  6.  However,  the  estimates  were  able  to  reproduce 
the  true  submarine  trajectory  reasonably  well. 
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An  area  that  remains  to  be  investigated  is  the  true  nature  of  the  errors  obtained  in  the  dynamic 
variables  x  ^  from  the  various  instrumentation  methods  to  be  discussed  in  the  following  section;  and 
how  this  affects  the  algorithm  as  constituted.  The  immediate  area  of  concern  is  the  cross  coupling 
that  will  be  obtained  between  the  dynamic  variables  due  to  the  non-orthogonality  of  the  instruments 
or  their  mountings.  Some  comfort  is  gained  in  that  the  white  noise  level  (Qj)  can  always  be  increased 
in  the  Kalman  filter  such  that  unmodeled  errors  reside  well  within  this  level.  This  generally  assures  a 
stable  computation  with  some  loss  of  optimality. 
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X.  NAVIGATION  AND  MOTION  MEASUREMENT  INSTRUMENTATION 


The  submarine  parameter  identification  program  is  dependent  on  motion  measurement  data 
of  high  resolution  over  several  minutes  of  maneuvering  flight  time.  Medium  to  high  accuracy  data 
are  required  within  a  restricted  frequency  spectrum  of  the  submarine  motion  during  the  maneuvers. 

The  complexity  of  the  on-board  motion  measurement  instrumentation  is  partly  a  function  of  the 
availability  and  accuracy  of  various  navigation  equipments.  Several  options  are  outlined,  including 
an  inertial  platform  navigator,  and  an  inertial  strap-down  navigator. 

A  possible  solution  to  the  instrumentation  problem  is  a  low  cost  inertial  navigator  with  a  verti¬ 
cal  channel  slaved  to  pressure  depth.  The  navigator  North  and  East  Schuler  axes  are  damped  by  the 
ship’s  EM  log  resolved  by  pitch  and  azimuth  through  a  MK  19  gyro  compass.  The  azimuth  gyro  of  the 
navigator  is  aligned  and  drift  corrected  by  feeding  back  a  comparison  of  platform  azimuth  with  the  MK 
19  azimuth. 

If  the  submarine  is  traveling  at  reasonably  constant  depth,  heading,  and  speed  for  approximately 
one  hour  the  inertial  navigator  can  become  aligned,  leveled,  and  referenced  to  the  average  local  ocean 
currents.  Before  starting  maneuvers  the  inertial  navigator  is  uncoupled  from  the  EM  log  and  the  gyro 
compass.  The  navigator  has  less  dynamic  error  than  the  other  equipment,  for  accelerating  maneuvers 
over  short  periods  of  time. 

The  submarine  performs  the  maneuvers  based  on  pseudo-random  pulse  sequences  applied  to 
the  bow  plane,  stern  plane  and  rudder.  Control  biases  are  added  to  the  random  pulse  sequence  to  avoid 
excessive  pitch  angles  and  depths.  The  control  surface  deflections  must  be  very  accurately  measured 
but  need  not  necessarily  be  accurately  controlled.  Manual  inputs  obtained  by  following  a  panel  meter 
display  are  satisfactory. 

After  several  minutes  of  maneuvering,  the  submarine  returns  to  its  original  depth,  heading,  and 
speed.  At  this  time  the  velocity  outputs  of  the  navigator  are  recorded  and  compared  with  the  EM  log 
and  gyro  compass  to  obtain  an  estimate  of  the  Schuler  error  build  up. 

Assuming  the  ocean  currents  are  constant,  the  inertial  platform  measurements  during  the 
maneuvers  are  with  respect  to  the  water  mass.  This  is  desirable  because  the  basic  equations  are  hydro- 
dynamic. 

The  platform  errors  are  mainly  at  the  Schuler  frequency  (a  period  of  84  minutes)  of  about 
0.001 25  rad/sec.  The  recommended  pseudo  pulse  sequence  has  a  minimum  frequency  of  about  0.02 
rad/sec  and  a  maximum  significant  frequency  of  about  0.6  rad/sec.  The  Schuler  errors  of  the  platform 
can  be  separated  in  the  frequency  domain  from  the  measured  motion  of  the  submarine  during  the  pseudo¬ 
random  maneuver.  This  insures  that  the  Schuler  cycle  can  be  estimated  and  corrected  in  order  to  remove 
ramps  and  bias  caused  by  the  Schuler  oscillation  during  the  maneuver. 

A  significant  instrumentation  problem  involves  converting  the  navigator  outputs  in  terms  of 
earth  fixed  coordinates.  North,  East  and  vertical  velocity  and  roll,  pitch  and  yaw,  to  body  coordinates 
of  u,  v,  w,  p,  q,  r  and  derivatives  thereof.  The  method  of  accomplishing  this  is  shown  in  Figure  3.  The 
coordinate  transformation  represented  by  Figure  3  must  be  accomplished  with  an  inaccuracy  of  less 
than  0. 1  %  preferably  0.05%. 


The  derivatives  can  be  obtained  by  a  combination  of  processing  and  smoothing  of  the  basic 
quantities  in  such  a  way  that  they  are  in  juxtaposition  with  each  other,  i.e.,  lagged  equivalently  in  time 
or  processed  by  equivalent  spectral  windows,  such  that  the  smoothed  quantity  is  the  integral  of  the 
smoothed  derivative. 


A  possible  alternative  to  the  standard  platform  type  inertial  navigator  is  the  strap  down  iner¬ 
tial  navigator.  The  strap  down  instrumentation  is  attractive  because  the  basic  measurements  are  al¬ 
ready  in  the  body  axes  coordinate  system.  The  requirements  of  coordinate  transformations  through 
roll,  pitch  and  yaw  still  exist  however,  since  the  strap  down  system  must  be  pseudo  leveled  with  respect 
to  the  gravity  vector.  Figure  4  shows  how  a  typical  vertical  gyro  may  be  simulated  by  body  mounted 
accelerometers  and  rate  gyros.  In  order  to  understand  Figure  4  the  following  equations  are  helpful. 


1) 


2) 


3) 


4) 


5) 


6) 


p  i 

q  =  R  8 

r  i 

0  P 

6  =  R" 1  q 

i  r 

1  0  -sin  8 

R  =  0  cos  0  cos  8  sin  0 

0  -sin  0  cos  8  cos  0 


1  sin  0  tan'  6  cos  0  tan  6 

RT1 

= 

0  cos  0  -sin  0 

0  -sin  0  sec  8  cos  <p  sec  8 

ax 

= 

•8x+au 

ay 

= 

gy  +  av 

az 

= 

gz  +  aw 

au 

= 

u  +  qw  -  rv 

av 

= 

v  +  ru  -  pw 

aw 

= 

w  +  pv  -  qu 
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The  vertical  gyro  simulation  as  shown  is  capable  of  a  nominal  roll  and  pitch  measurement  simil¬ 
ar  to  that  of  an  actual  gyro.  It  is  reasonably  accurate  under  normal  cruise  conditions.  The  gyro  leveling 
loops  closed  by  gain  Kp  and  Kr  may  be  momentarily  removed  in  which  case  the  gyro  simulated  is  a  free 
gyro.  The  roll  axis  gain  Kr  is  typically  removed  during  turns.  If  the  rate  gyro  and  accelerometer  outputs 
are  encoded  for  transmission  into  a  digital  computer,  the  mechanization  shown  becomes  strictly  computa¬ 
tional.  A  simulated  free  vertical  gyro  during  maneuvers  may  be  corrected  for  earth  rate  and  motion  over 
the  earth.  In  this  case  the  simulation  begins  to  resemble  a  strap  down  inertial  system. 

The  next  step  in  strap-down  motion  measurement  instrumentation  is  to  indicate  the  body  axis 
velocity  computations.  Figure  5  shows  a  method  of  instrumenting  u,  v,  and  w  where  u  is  continually 
compared  with  the  EM  log.  The  u,  v,  and  w  components  are  then  resolved  through  roll  and  pitch,  inte¬ 
grated  and  compared  with  depth.  Errors  between  computed  depth  and  actual  depth  are  used  to  correct 
u,  v  ,  and  w  and  the  gravity  computation.  Other  outputs  such  as  and  VgaS|.  can  be  compared 

with  SINS.  and  V£ast  can  be  integrated  in  terms  of  Latitude  and  Longitude  and  these  can  be 

compared  with  outputs  from  Omega  or  SINS  etc.  The  total  mechanization  of  Figures  4  and  5  can  be 
computational.  The  gains  Kp,  K,.,  K) ,  K2,  K3,  and  K4  can  be  based  on  a  covariance  computation  as 
per  the  optimal  Kalman  filter. 

p,  q,  r  and  p,  q,  r  can  be  obtained  with  good  accuracy  by  using  a  rate  integrating  gyro  with  feed¬ 
back  between  pickoff  and  torquer.  Inaccuracies  occur  due  to  scale  factor,  misalignment,  drift  of  the 
non-g,  g  and  g2  type,  output  axis  coupling,  earth  rate,  motion  over  the  earth,  and  coneing.  Nearly  all 
these  effects  can  be  compensated  by  careful  measurement  of  the  appropriate  gyro  characteristics  be¬ 
fore  installation  and  mechanizing  the  necessary  computational  feedbacks  between  the  gyros, acceler¬ 
ometers  and  other  motion  measuring  equipment.  Of  course  the  gyro  parameters  are  random  and  vary 
a  certain  amount  from  minute  to  minute,  day  to  day,  or  turn-on  to  turn-on,  consequently  the  compensa¬ 
tion  is  never  perfect.  Lack  of  compensation,  mount  vibration,  bearing  rumble  and  rotor  RPM  varia¬ 
tions  can  be  considered  as  gyro  noise.  Figure  6  avoids  these  complications  but  shows  the  basic  p  and 
p  mechanization.  The  computation  can  be  made  accurate  from  0  to  50  Hz  (with  corrections).  Addi¬ 
tional  filtering  (not  shown)  to  eliminate  ship’s  vibration  and  noise  is  desirable  to  prevent  aliasing  of 
the  data  before  any  digital  processing  is  attempted.  The  same  filters,  if  required,  should  be  used  on  the 
accelerometers.  The  filters  should  have  negligible  phase  shift  or  attenuation  from  zero  to  0.5  Hertz. 

(This  may  not  be  adequate  for  model  submarines). 
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pmeas  pmeas 

Figure  6.  Basic  p  and  p  Measurement 
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XI.  NOISE  FREE  PARAMETER  ESTIMATE  TRAJECTORIES 


This  section  contains  typical  trajectories  of  the  estimates  of  the  parameters,  or  coefficients, 
under  the  data  processing  algorithms  described  in  Section  VI.  The  parameter  trajectories  shown  in 
this  section  are  computed  with  the  assumption  of  sensor  bias,  but  no  sensor  noise.  The  sequence  is: 

Axial  Force  Parameters 
Rolling  Moment  Parameters 
Pitching  Moment  Parameters 
Normal  Force  Parameters 
Lateral  Force  Parameters 
Yawing  Moment  Parameters 

All  trajectories  of  the  estimates  have  been  normalized  so  that  correct  estimation  yields  an 
ordinate  of  unity. 

A  common  format  is  used  throughout;  the  entire  set  of  data  is  processed  twice.  On  the  first 
pass,  five  loops  are  conducted  through  each  equation  from  time  equal  to  zero  and  300  seconds.  On 
the  second  pass  only  one  loop  is  conducted.  This  corresponds  to  a  maximum  NP  of  2  where  NL 
equals  5  for  NP  =  1  and  NL  equals  1  for  NP  =  2  in  Figure  1 .  The  plots  for  the  hydrodynamic  co¬ 
efficients  show  the  entire  estimation  history  in  the  algorithm  of  Figure  1 . 

Since  estimates  for  the  bias  errors  (in  the  dynamic  state  variables)  are  obtained  from  all  the 
loops  of  every  equation,  a  different  format  is  employed  in  the  plots  of  these  quantities.  Only  the 
estimation  response  from  the  final  loop  for  a  given  NP  is  shown  in  these  figures. 

As  the  digital  plot  machine  does  not  provide  subscripts,  it  is  necessary  to  interpret: 

XQQ  =  Xqq,  etc. 

It  should  be  emphasized  that  the  notation  used  in  the  plots  are,  in  some  cases,  not  the  same  as 
Tables  1  through  6.  The  subscripting  of  the  hydrodynamic  parameters  in  the  tables  is  that  of  Reference  3, 
which  is,  for  some  parameters,  overly  abbreviated.  The  actual  nature  of  the  hydrodynamic  effect  is  there¬ 
fore  obscured.  For  this  reason  the  notation  used  in  the  plot  figures  has  been  modified  so  that  the  true 
nature  of  the  hydrodynamic  effect  is  indicated.  The  presentation  of  the  parameters  in  the  plots  is  in 
the  same  order  as  that  of  the  tables. 

The  phenomena  of  rapid  fluctuation  observed  in  the  yaw  curves  near  the  beginning  of  pass  2  is 
not  readily  understood,  since  one  would  not  expect  it  to  occur.  It  should  be  pointed  out  that  on  pass  2 
the  covariance  matrix  for  the  parameters  is  restarted  and  the  measurement  noise  level  is  significantly 
lowered.  Therefore  the  Kalman  gains  are  “wide-open”  to  allow  large  deviations.  Note  also  that  the  roll 
and  pitch  angular  acceleration  rate  errors  are  not  estimated  at  the  time  that  the  estimates  of  the  param¬ 
eters  dip.  These  may  account  for  the  dip  that  the  estimates  show  in  the  interval  70<t<  150  seconds. 
Consider  also  that  unique  determination  of  the  yaw  coefficients  is  heavily  linked  to  blowing  ballast,  which 
occurs  at  150  seconds  when  the  rapid  convergence  occurs.  This  effect  accounts  for  the  convergence 
to  correct  values  observed  at  t>  150  seconds.  The  covariance  matrix  for  the  parameters  did  show  a 
rapid  decrease,  so  that  a  rapid  fluctuation  of  the  estimates  is  not  unexpected. 
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XII.  PARAMETER  ESTIMATE  TRAJECTORIES  WITH  NOISE 


This  section  presents  trajectories  showing  the  convergence  of  the  normalized  estimates  of 
the  various  parameters,  or  coefficients,  when  the  sensor  data  were  corrupted  by  bias  plus  white 
noise. 

The  comments  of  Section  XI  on  the  sequencing  of  the  equations  and  interpretation  of  the 
subscripts  apply  equally  to  this  section.  As  in  Section  XI,  the  dynamic  state  variables  are  also 
shown  in  a  normalized  form. 
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APPENDIX  A 

DISCUSSION  OF  THE  COMPUTATIONAL  METHOD 


Before  proceeding  it  is  necessary  to  point  out  that  the  symbols  used  in  this  section  are 
for  the  purpose  of  discussion,  and  are  not  to  be  confused  with  the  use  of  similar  symbols  in  the 
body  of  the  report. 

The  general  algorithm  presented  in  Section  VI  is  difficult  to  explain  in  terms  of  a  simple 
system  —  since  it  was  specifically  developed  to  handle  complex  systems  that  are  structured  with 
a  large  number  of  variables  (exceeding  100).  However,  a  conceptual  basis  that  will  facilitate 
understanding  of  the  general  algorithm  can  be  obtained  via  a  simple  example. 

Let  z  be  a  scalar  variable  whose  time  derivative  is  governed  by  the  following  equation, 

z  =  gi  (z,  uc)  +  Cjg2  (z,  uc)  +  c2g3  (z,  uc) 

where  uc  is  a  control  variable  and  Cj ,  c2  are  unknown  coefficients;  gj  (z,  uc),  i  =  1,2,3  denotes  a 
general  function  that  is  dependent  upon  z  and  uc.  It  is  assumed  that  z  and  z  are  directly  instru¬ 
mented  or  constructed  from  measurements  that  are  a  function  of  them.  It  is  necessary  to  recog¬ 
nize  that  z  and  z  are  not  perfectly  known. 

In  the  following  the  central  idea  is  that  the  error  characteristics  of  the  estimates  for  z  and 
Let  z,  z  denote  the  estimated  values  of  z  and  z;  then  the  errors  in  z  and  z  are  de- 


-  z 

-  z* 

the  nature  of  z”and  z  is  usually  given  in  terms  of  their  time  derivatives.  To  avoid  confusion,  let 
Xj  =Tand  x2  =  z.  If  the  trivial  case  where  Xj  and  x2  are  constant  is  assumed,  then 

Xj  =0 

x2  =  0. 

Note  that  in  this  situation  z"#  z  if  for  example  z  and  z  are  measured  separately  by  different  in¬ 
struments.  These  considerations  also  apply  to  the  parameters  cr  and  c2 .  From  an  a  priori  under¬ 
standing  of  the  physical  situation  represented  by  the  differential  equation  for  z,  the  nature  of  the 
uncertainty  in  Cj  and  c2  can  also  be  ascertained. 

If  a  priori  values  for  Cj  and  c2  can  not  be  guessed,  the  choice  of  zero  becomes  appropriate. 
As  above,  the  estimated  values  for  Cj  and  c2  are  denoted  by  and  £2 .  If  it  is  known  that  Cj  and 
c2  are  constant,  then  the  errors  in  the  estimates  are  characterized  as  follows: 

^  =  0, 

c"2  =  0. 


z  be  specified, 
noted  by 


z  = 


z  = 
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Let  Cj  =  x3  and  ~c2  =  x4)  then  the  error  dynamics  for  xj,  i  =  1,  2,  3,  4  can  be  summarized  in 
vector  form 


In  the  following  x  will  be  referred  to  as  the  state  vector.  A  more  general  linear  stochastic  model 
for  x  may  of  course  be  possible, 

x  =  Fx  +  uw 

where  F  is  a  matrix  and  uw  is  a  white  noise  vector.  This  is  the  familiar  form  for  error  dynamics 
in  which  Kalman  filtering  is  applied.  For  the  example  above  F  and  uw  are  a  zero  matrix,  and  zero 
vector  respectively. 

Let  the  differential  equation  for  z  be  written  in  terms  of  the  estimated  quantities, 

A  A  A  .  A  /A  A.  A  .A  A. 

Z  =  gi  (z,  uc)  +  C!g2  (z,  Uc)+c2g3  (z,  Uc). 

Note  that  A  appears  over  uc  since  it  also  is  an  imperfectly  measured  quantity.  For  the  following 
it  will  be  assumed  that  uc  will  have  only  a  white  noise  component  of  error.  It  is  readily  recognized 
that  if  the  estimated  quantities  are  not  perfect  then  the  equation  above  will  not  be  satisfied.  An 
error  signal  can  now  be  defined, 

~  _  A  A  A 

y  -  z  -  gi  -C!g2  -c2g3. 

For  ease  of  presentation  the  following  definition  has  been  introduced 

gi  =  gi(£,uc). 


This  equation  represents  exactly  how  the  actual  error  signal  would  be  computed.  The  problem  now 
is  to  express  y- in  terms  of  the  basic  error  vector  x.  A  linear  approximation  is  of  course  desired  so 
that  the  Kalman  filter  algorithm  can  be  employed  to  computed  gains  that  multiply  y-  and  obtain  a 
correction  for  the  estimated  state  vector  x.  Since  z  =  z  +  'z,  y"  can  be  expanded  about  the  true 
variables, 


Z  '  gl  (Z>  Uc)  -  Cjg2  (z,  uc)  -  c2g3  (z,  uc) 

.t-  9gi  _  9gj  _  _  dg2  - 

+  z '  —  z  -  —  uc  -  clg2  -Cl  — 
9z  9uc  w  9z 


3g2  ~ 

z  -  c,  -  u. 

19uc  ( 


‘  c2  g3  '  c2  —  '  z  -  c2  —  ur  +  higher  order  terms. 

9z  9uc  L 

By  hypothesis  the  first  four  terms  above  add  up  to  zero;  after  collecting  terms  y”  can  be  put  into 
the  following  form, 
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y  =  Mx  +  v  +  higher  order  terms 
where  M  is  a  single  row  matrix  whose  elements  are  given  by 


mi 

dgi 

3z 

m2 

=  1 

m3 

=  -g2 

m4 

=  -g3 

V 

=  JiSL 

[3uc 

dg2  c  dg3 


Since  the  evaluation  of  g^  and  the  associated  partial  derivatives  can  only  be  done  with  estimated 
quantities  an  additional,  higher  order,  error  is  introduced  into  the  linear  expansion  fory”.  Use  of 
the  Kalman  filter  algorithm  implicitly  ignores  these  higher  order  terms.  However,  the  effect  of 
these  terms  can  be  accounted  for  by  computing  the  Kalman  gains  with  a  larger  level  of  measure¬ 
ment  noise. 


Since  the  partial  derivatives  and  M  elements  (m^)  are  functions  of  the  states  which  are  to  be 
estimated,  multiple  computation  of  the  linearized  Kalman  filter  will  be  required.  The  phrase  “com¬ 
putation  of  the  Kalman  filter”  implies  the  computation  of  the  error  covariance  matrix,  optimal 
gains  and  estimation  of  the  states  for  the  time  interval  that  data  are  available.  Throughout  the 
remainder  of  this  report  this  computational  sequence  will  be  designated  as  a  sweep  or  loop. 
Typically  Cj  and  c2  would  be  initialized  at  zero,  so  that  the  first  time  the  filter  is  iterated  (looped) 
rrij  and  the  computation  of  the  noise  variance,  E  |v2}  ',  would  be  incorrect.  It  is  stipulated  that  a 
loop  will  be  made  with  one  set  of  values  for  the  q’s.  This  is  an  important  feature  of  the  algorithm; 
generally  it  seems  natural  that  as  soon  as  improved  estimates  of  Cj  are  obtained  the  more  accurate 
computation  of  nij  should  be  made.  This  notion  is  incorrect  because  if  different  values  of  Cj  are 
used  in  the  computation  of  m^,  for  any  sweep,  a  high  degree  of  false  observability  is  inserted  into 
M  and  correspondingly  into  the  error  covariance  matrix.  The  covariance  matrix  becomes  unre¬ 
liable  and  divergence  is  likely  to  occur.  The  solution  taken  here  is  to  loop  the  Kalman  filter  with 
different  levels  of  measurement  noise  -  starting  at  some  high  level  and  decreasing  it,  possibly  in 
some  linear  or  quadratic  manner.  The  q  are  of  course  held  constant  for  the  sweep.  The  corrections 
for  obtained  from  the  loop  are  stored  into  a  separate  array  and  used  to  update  the  Cj  array  at  the 
completion  of  the  current  sweep. 

This  technique  can  be  readily  extended  to  the  multidimensional-nonlinear  parameter  esti¬ 
mation  problem.  The  procedure  is  the  same  as  that  followed  here;  only  in  this  case  M  will  be  a 
matrix  with  the  number  of  rows  equal  to  the  number  of  differential  equations  for  the  dynamical 
system. 


The  method  presented  here  is  similar  to  the  traditional  extended  Kalman  filter  applications 
to  estimation  of  nonlinear  systems,  in  the  sense  that  the  error  covariance  of  the  dynamic  variables 
(as  distinguished  from  system  parameters)  is  part  of  the  identification  procedure.  The  notion 
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introduced  here  is  that  the  nature  and  accuracy  of  the  actual  measurements  made  on  the  system 
uniquely  establish  the  error  covariance  matrix  for  the  dynamic  variables.  Thus,  a  velocity  meter 
and  accelerometer,  for  all  practical  purposes,  alone  establish  the  accuracy  (biases  and  covariance) 
of  velocity  and  acceleration,  and  little  or  nothing  is  gained  by  incorporating  the  linearized  differen¬ 
tial  equation  for  the  error  behavior  of  the  nonlinear  system,  as  is  done  in  standard  extended  Kal¬ 
man  filtering.  It  is  essentially  hypothesized  that  introduction  of  the  equations  of  motion  as 
measurements  on  the  error  space  will  be  adequate  to  obtain  a  stable  unbiased  estimation  algorithm. 

There  might  be  a  tendency  to  call  this  method  “Extended  Recursive  Least  Squares  Esti¬ 
mation”;  it  is  felt  that  the  term  “Least  Squares”  is  misleading  since  traditionally  a  direct  matrix 
inversion  (nonsingular)  is  associated  with  these  “Least  Squares”  methods.  The  Kalman  filter  is  an  ortho¬ 
gonal  projection  of  the  measurements  into  the  state  (error-parameter)  space  which,  taking  all 
error  characteristics  into  account,  yields  unbiased  estimates,  and  a  reliable  covariance  matrix  which 
indicates  what  parameters  are  poorly  estimated  and  to  what  extent  they  are  independent.  It  should 
also  be  mentioned  that  if  an  algorithm  is  unbiased,  it  minimizes  a  large  number  of  error  criteria 
(of  which  least-square  is  only  one). 

It  is  well  known  that  the  traditional  “Least  Squares  Method”,  which  ignores  white  noise 
and  bias  errors,  gives  poor  results  (4,  5,  6).  The  work  of  Kalman  (2)  has  elevated  the  level  of 
thought  towards  stochastic  filtering,  and  it  is  not  appropriate  to  refer  to  any  applications  of  his 
ideas  in  terms  of  ad  hoc  principles  such  as  least  square  or  maximum  likelihood. 
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